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We compare two calculations of the particle density in the superfluid phase of the classical XY 
model with a chemical potential p in 1+1 dimensions. The first relies on exact blocking formulas 
from the Tensor Renormalization Group (TRG) formulation of the transfer matrix. The second is 
a worm algorithm. We show that the particle number distributions obtained with the two methods 
agree well. We use the TRG method to calculate the thermal entropy and the entanglement entropy. 
We describe the particle density, the two entropies and the topology of the world lines as we increase 
p to go across the superfluid phase between the first two Mott insulating phases. For a sufficiently 
large temporal size, this process reveals an interesting fine structure: the average particle number 
and the winding number of most of the world lines in the Euclidean time direction increase by one 
unit at a time. At each step, the thermal entropy develops a peak and the entanglement entropy 
increases until we reach half-filling and then decreases in a way that approximately mirror the ascent. 
This suggests an approximate fermionic picture. 

PACS numbers: 05.10.Cc,05.50.+q,11.10.Hi,64.60.De,75.10.Hk 


I. INTRODUCTION 


The classical XY model in one space and one Euclidean 
time (1+1) dimensions plays an important role in 
the field theoretical approach of condensed matter 
phenomena and is prominently featured in standard 
textbookJil®. This model provides the simplest example 
of a Berezinski-Kosterlitz-Thouless transitiorf^ and it 
may be used as an effective theory for the Bose- 
Hubbard modeP. More generally, its associated quantum 
Hamiltonian of Abelian rotors appears in many different 
contexts such as the formulation of Abelian lattice gauge 
theorie^, Jose phson junctions array^ and cold atom 
simulat or 

When a chemical potential p is introduced, the model 
displays a rich phase diagram depicted in Ref. [101 
For sufficiently small /3, the inverse temperature of the 
classical model, if we increase p, we go from a Mott 
insulating (MI) phase where the average particle number 
p remains zero until it reaches a superfluid (SF) phase 
where p starts increasing with p. This proceeds until 
p reaches one per site and we enter in a new MI phase 
where it stabilizes at this value. For /? small enough, this 
alternation of MI and SF phases repeats several times. 
In suitable coordinated^, the phase diagram is similar 
to what is foun d for the one-dimensional quantum Bose- 
Hubbard modefi^li^l. 


A simple way to locate approximately the SF phase 
consists in calculating the thermal entropy, defined 
precisely in Eq. (161, on a lattice with a large enough 
temporal size. An example is shown in Fig. The 
central line at /3=0.1, where we do calculations hereafter, 
covers the first SF phase and the MI phases with p = Q 
and 1. As we increase p the thermal entropy goes 


through a sequence of peaks culminating around In 2, 
signaling level crossings that will be explained below. 
The number of peaks equals the number of sites in the 
spatial direction. This is clearly a finite size feature, while 
the notion of phase used above should be understood in 
the limit of an infinite number of sites. 
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FIG. 1. (Golor online). Intensity plot for the thermal entropy 
of the classical XY model on a 4 x 128 lattice in the fi-p plane. 
The dark (blue) regions are close to zero and the light (yellow 
ochre) regions peak near In 2. The MI phase with p=0 is below 
the lowest light band, the MI phase with p—1 is above the 
highest light band and there is a single SF phase in between 
the two MI phases. 


In this article, we describe microscopically the rich 
sequence of changes which occurs as we move across the 
SF phase along a line of constant /? as described above 
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and illustrated in Fig. The figure can be embedded 
in the phase diagrams of Ref. HOl the lower and upper 
light lines ending at the tips of the the two MI phases. 
The fine structure that we report here is first studied for 
a small spatial size of four sites and then for larger sizes. 
At infinite spatial size, there are three phases in Fig. [l] 
the SF phase approximately covering the four light bands 
with larger values of the thermal entropy and the three 
darker regions in between. We also report about the 
numerical methods that we developed in this process. 

The specific classical XY model used in this article 
is a planar version of the Ising model with a 0(2) 
symmetry, on a x Lt lattice. The notations used 
later to characterize the model and the basic numerical 
methods are provided in Sec. [IT] Most of the calculations 
done in this article will rely on the tensor renormalization 
group (TRG) methocP^ which can be used to write exact 
blocking formulas for the model considerecP^®^ . In Sec. 
|II B[ we remind how the method can be applied to the 
calculation of the transfer matriisP^. 


It should be noticed that the presence of fi causes the 
action to be complex. This sign problem prevents the use 
of Monte Carlo simulations when /r is too large to rely on 
reweighing methods. However, using a Fourier expansion 
of the Boltzmann weightJi^, it is possible to find a 
formulation of the partition function in terms of world 
lines with a positive weight as long as /i is real. This 
positivity allows statistical samplin^i^ using a classical 
version of the worm algorithirPHl. The computational 
methods used for the worm algorithm are briefly reviewed 
in Sec. Ill Cl 

The TRC method can be used for arbitrary complex 
value^^ of P and /r. The only sources of error are the 
truncations of the infinite sum to finite ones as required 
for the numerical treatment. These truncations take 
place in the original formulation of the partition function 
and also at each step of the coarse-graining process. A 
first check of the agreement between the TRC and the 
worm methods is provided in Sec. [m] where we compare 
particle number density histograms and show that they 
agree well. 

In the SF phase, the model is gapless in the the limit 
Lx —> oo. At finite Lx and small /3, the gap is expected 
to scale like l/T^- This follows from a non-relativistic 
dispersion relation and can be justified using degenerate 
perturbation theory. An invaluable method to study 
the long range correlations in near gapless situations 
is to calculate the entanglement entropy. Interestingly, 
it seems possible to measure the entanglement entropy 
in many-body systems implemented in cold atoms^SHSSl 
This includes studies of rather small Lx systems. The 
entanglement entropy needs to be distingui she d from the 
thermal entropjEH. For this reason, in Sec. IV we discuss 
these two entropies using the transfer matrix formalism. 
We u se a l attice version of the setup of Calabrese and 
CardjESHni for 1+1 dimensional nonlinear sigma models. 
We consider the case of finite, but often large, Lj. A 
finite Lt introduces a temperature proportional to 1/At, 


distinct from 1//3 used in the classical formulation, and 
therefore we have a thermal density matrix. In this 
context, the relation between the two entropies is a rather 
open topic of investigatiorP^. Numerical calculations of 
the related Renyi entropy of the classical XY model, 
without chemical potential, were presented in Ref. 1321 

With the TRC method, we approximate the reduced 
density matrix by a finite dimensional matrix which can 
be diagonalized numerically and we do not need to use 
the replica trick as in Refs. [29l[30l and[32] We then show 
how to use the TRC method to express the entanglement 
entropy for a bipartition of the system in two subsystems. 
We show that for small Lt, the thermal entropy is larger 
than the entanglement entropy but the thermal entropy 
becomes larger as Lt is increased. 

In Sec. |Vj we use degenerate perturbation theory to 
get an approximate idea of the large structure (location 
of SF phases) and fine structure (changes of p and 
entanglement entropy across one SF phase) of the phase 
diagram. We relate in some approximate way, the 
eigenvectors of the transfer matrix with particle number 
n to world-line configurations with a winding number 
which is also n. We then calculate numerically the 
average particle number density, the thermal entropy and 
the entanglement entropy for values of p spanning the 
first SF phase as explained above. 

For sufficiently large Lt, the results show an interesting 
fine structure: the particle number and the winding 
number of (most of) the world lines increase by one unit 
at a time as we keep increasing p. At each step, the 
thermal entropy develops a peak and the entanglement 
entropy increases with p until p reaches half-filling. As 
we keep increasing p beyond this value, the entanglement 
entropy decreases in a way that approximately mirrors its 
ascent. This approximate symmetry can be justified by 
noticing that if we interchange the occupied links with 
the unoccupied links in the world lines, we transform 
a configuration with particle number n to one with a 
particle number Lx — n. This approximate particle- 
hole transformation can be reformulated in the context 
of degenerate perturbation theory and suggests that a 
fermionic description is possible. 

Our results are summarized in Sec. IVII where we also 
briefly discuss work in progress. We suggest ways to 
reduce the small truncation errors reported in Sec. |III| 
and to interpret the approximate particle-hole symmetry 
found in Sec. |V] We also briefly discuss t he re lationship 
of our work with Polyakov’s loop studie^^^*^ and with 
recently proposed cold atom experiments^. 


II. NOTATIONS AND NUMERICAL METHODS 

In this section, we introduce the classical XY model 
with a chemical potential. We describe the two numerical 
methods used in this article. The first is the TRC 
which relies on coarse-graining, the second is a worm 
algorithm which relies on sampling. This section contains 
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important concepts used in the next sections, such as 
the representation of the partition function in terms of 
integers labeling Fourier modes attached to the links, also 
called bounds, of the lattice, their use in the transfer 
matrix formulation and their graphical illustration as 
world lines. The section also contains technical details 
about the numerical implementations that are given for 
completeness. 


A. The model 


We consider the classical XY model, sometimes called 
the 0(2) model, with one space and one Euclidean time 
direction, and a chemical potential /i. One can interpret 
/j, as the imaginary part of a constant gauge field in the 
temporal direction. The sites of the rectangular lattice 
are labelled (x,t) and the unit vectors denoted x and 
t. The line segments joining two nearest neighbor sites 
are called links. The total number of sites is x Lt- 
In the following, we are typically interested in the case 
Lt >> Lx- We assume periodic boundary conditions in 
space and time. The partition function reads 


Z = 


n dO{x,t) s 

27T 

(x,t) 


( 1 ) 


with 

{x,t) 

y ^ C0s(d(a;_|_l ( 2 ) 

{x,t) 


In all the numerical calculations done in this article, we 
consider space-time isotropic couplings j5x = Pt = P- 
However, if we set Px to zero, the model becomes a 
collection of decoupled solvable models. For analytical 
purposes, when P is small, it is sometimes convenient to 
first consider the solvable case Px=0 and then restore the 
isotropic situation Px = Pi = P perturbatively (see Sec. 


As explained in Refs, [ini HU and m one can use the 
Fourier expansion of in terms of modified Bessel 

function of the first kind and then integrate out the 
d{x,t) variables. The Fourier indices associated with the 
links coming out of the site {x^t) in the space and time 
directions are denoted n(^x,t),x and i respectively. 
They can be interpreted as currents or particle numbers 
passing through the links. The partition function can 
then be expressed as a sum of product of Bessel functions: 


{n} (x,t) 

The Kronecker delta function in Eq. @ ensures the 
local current conservation and the terms in the partition 
function can be interpreted as current loops (also called 


world lines) which can then be statistically samplecP^. 
For a system with periodic boundary condition in both 
space and time directions, as considered in this paper, 
the world line can wind around in both directions. 
The winding numbers are important to understand the 
superfluid properties of the systenP^. 

)S=0.1; )U=3.0 



X 

FIG. 2. (Color online). Graphical representation of an 
allowed configuration of {n} for a 4 by 32 lattice. The 
uncovered links on the grid have n=0, the more pronounced 
dark lines have |n| = l and the wider lines have n=2. The 
signs of the links with |n|=l lines are discussed in the text. 
The large dots on the horizontal boundaries (red) need to be 
identified in pairs with the same x coordinates. Similarly, the 
slightly smaller dots (blue) on the vertical boundaries have to 
be identified in pairs with the same t coordinate. 


An typical allowed configuration for /3=0.1 and /r = 
3 is shown in Fig. for Lx=^ and L(=32. Sites at 
the boundary should be identified as explained in the 
figure caption. The uncovered links on the grid have 
n=0, the more pronounced dark lines have |n| =1 and 
the wider lines have n=2. A discussion regarding the 
sign convention and the spatial winding number of Fig. 
[^can be found in Appendix [A| 

This configuration can be used to visualize a transfer 
matrix that connects consecutive time slices. For 
instance, in Fig. the time slice 5 represents a transition 
between |1100) and |0200) and its relative statistical 
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weight can be obtained from a transfer matrix that 
will be discussed in Sec. |IIB The configuration was 
generated using a sampling method designed in Ref. ITOl 
and briefly discussed in Sec. m 


transfer matrix elements are zero unless the sum of the 
two sets of indices (respectively denoted n and n' in Eq. 
(§) are equal. In other words, the transfer matrix is 
block diagonal in each particle number sector 


B. TRG approach of the transfer matrix 


Lx Lx 

n = '^n^ = ■ (7) 

X — 1 X — 1 


As explained in Ref. (TUI the partition function can be 
expressed in terms of a transfer matrix: 

Z = TrT-^* . (4) 


The matrix elements of T can be expressed as a product 
of tensors associated with the sites of a time slice (fixed t) 
and traced over the space indices. To make the equations 
easier to read, we use the notations for the time indices 
in the past f), the primed symbol for the time 

indices in the future and fix for the space indices 

i'n[x,t),x)- The matrix elements of T have the explicit 
form 


with 


T 




E rp{l,t) rp'' 


d2.i) 

'nih2n2n2. 


nin2...nL^ 




( 5 ) 


fix-ifixrixn'^ 


\J{PiVk if^i) expiniux + n'x)) 

1 iPx)Ih^ i/3x)Sn^-i+n^,n^+n'^ .( 6 ) 


In Eq. the trace is over the temporal indices, 

represented graphically by vertical lines in Eig. while 
the spatial indices (horizontal links in Fig. are summed 
over as described in Eq. ([^. The Kronecker delta 



/?! 


When the chemical potential \i is zero, there is a 
charge conjugation symmetry which allows us to change 
the sign of all the n’s in all the temporal sums without 
affecting the final results. Given the rapid decay of the 
modified Bessel function when the index n increases, 
good approximations can be obtained by replacing the 
infinite sums by sums restricted from —Umax to rimax- 
When the chemical potential is nonzero, it is more 
efficient to shift the range in the same direction as the 
sign of the chemical potential. In general, we call Dg the 
number of states kept after truncation. In the following, 
we will use a coarse-graining procedure for the transfer 
matrix where we repeatedly apply a truncation at the 
same value Dg. 

For numerical purposes, we reduce the si ze of the 
transfer matrix by using a blocking procedur^SHH j|- 
consists in iteratively replacing blocks of size two in the 
transfer matrix by a single site using a higher order 
singular value decomposition. During the truncation, 
for a given pair of sites, we replace the direct product 
X Dg matrix by a Dg x Dg matrix. The process is 
illustrated in Fig. The Y-shaped parts in the bottom 
part of the figure and their mirror versions in the top 
part schematically represent this truncation. The full 
technical details can be found in Refs. M and [ini 



FIG. 4. Graphical representation of the coarse graining 
truncation of the transfer matrix described in the text. 


FIG. 3. Graphical representation of 

qr 

function in Eq. ®) is the same as in Eq. ([^ and reflects 
the existence of a conserved current as discussed above, 
that we will call “particle number” in the following. For 
periodic (spatial trace) or open (zero spatial indices at 
both ends), the local conservation law implies that the 


C. The worm algorithm 

Because of the conservation law in Eq. ([^, the terms 
of the partition function can be inter pre ted as current 
loops which can be statistically sampled^. The sampling 
procedure for one complete worm algorithm update goes 
as follows. We pick randomly a site on the lattice, then 
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pick randomly a neighboring direction in the (positive or 
negative) spatial or time direction, change the current n 
to n ± 1 (depending on the sign of the change) and move 
to the neighboring site if the above update is accepted 
(the details of the accept-reject procedure are given in 
the Appendix A of Ref. SaD. We repeat the above 
procedure until we come back to the original random site. 
Because of the conservation law, a particle number can 
be attributed to a configuration. It can be calculated by 
summing the temporal n’s between any two time slices. 
The particle number distribution, can be calculated by 
generating a large number of configurations according to 
the above procedure. 

An easily controllable source of error is the limited 
statistics. There is also an unavoidable truncation error. 
Following Ref. (THl |n| is constrained to be not larger 
than 20 in Eq. However, given that /2o(0.1) exp(3 x 
20)//o(0.1) ~ 4 X 10“^®, this is negligible for our 
calculations. By construction, there is a single worm 
and the algorithm cannot generate disconnected loops. 
This is not believed to be a significant source of error. 
Comparisons at small volume where the truncation errors 
are controllable^ indicate that the worm algorithm is 
statistically exact. 


III. PARTICLE DENSITY CALCULATIONS 


In this section, we show how to calculate the average 
particle density using the TRG formulation. We then 
apply the method numerically and compare the results 
with the ones obtained with the worm algorithm. The 
particle number conservation can also be exploited in the 
TRG approach. For the initial one-site tensor we 

have 

^ = ^(no. + . (8) 


Here, Ux and are the particle number associated with 
the time indices of the original tensor with the 

tensor indices omitted. Consequently, we can associate 
a particle number n{i) with each eigenvalue of the 
transfer matrix T. The average particle number density 
is an extensive quantity defined as 


1 d\nZ 
ZxLt dji 


( 9 ) 


From the expression of Z in terms of the transfer matrix 
and the cyclicity of the trace, we have 


dZ/dii = LtTr(T'T^*-i), (10) 


where T' = dT/d^j, can be calculated by using the 
chain rule in Eq. ([^. This can be achieved iteratively 
by defining an “impurity” tensor initialized with the 
derivative of the initial tensor and then blocked and 
symmetrized with the original “pure” tensor. This 
guarantees the recursive replacement of by 


in the transfer matrix as prescribed by the chain rule. 
This procedure can be shortcut if we know the particle 
number n{i) associated with each eigenvalue as 
discussed above. We can then write 

IdlnZ E.A^n(z) 

Lt 5m ^ ^ 

In practice, finding the particle number associated with 
the eigenvalues is not completely straightforward. It 
requires to keep track of the particle number in the 
projected basis or to write the blocking algorithm sector 
by sector. The enforcement of the conservation law 
in the coarse-graining process makes the numerical 
calculation more stable. The tensor elements violating 
the conservation law are exactly zero. If we only handle 
non-zero elements, we can reach relatively larger values 
of the dimension Dg in the truncation procedure and get 
more accurate results. These improvements have been 
used in the numerical calculations done below. 

Knowing the n{i) associated with Ai also allows us to 
define a probability P{n) for the particle number n: 

P(.n)= Y. (12) 

i:n{i)—n i 

These probabilities can also be calculated directly from 
histograms obtained with the worm algorithm. Both 
methods can be used to calculate and compare the 
average particle number density using 

p=^YnP{n). (13) 

^ n 

We studied numerically the distribution of P{n) for 
various values of /i spanning a range covering the 
boundary of the MI phase and the SF phase. The 
other parameters are kept fixed at Dg = 201, /? = 0.1, 
Lx = 32, and Lt = 128. When p = 2.8 and 2.85, 
the distribution bears one-bin structure with n = 0, 
corresponding to the MI phase. When p = 2.9, 2.95 
and 3, the distribution carries more bins, corresponding 
to the SF phase, which shows that the phase transition 
occurs in the range p = (2.85, 2.9) as illustrated in Fig. 

We then compared the distributions obtained with the 
TRG and the worm algorithm for p=3, which is near the 
middle to the SF phase. The results are shown in Fig. 
1^ To the best of our knowledge, the errors associated 
with the worm calculations are purely statistical. The 
errors made with the TRG are due to the repeated 
truncation to Dg indices. We have used iAs=101, 201 
and 301 and the variations are comparable to the worm 
errors. Overall, the results from L>g = 301 are very 
close to the worm calculation. The particle number 
distribution calculated from the worm histograms are 
from averaging over configurations generated with twelve 
different initial random seeds. With each seed, four 
million equilibrium configurations were generated. With 
the spatial dimension increasing, the accuracy becomes 
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FIG. 5. (Color online). The particle number distribution 
P{n) with ^ taking different values at the boundary between 
the MI and SF phases. For /r = 2.80 and 2.85, we only have 
n = 0. For fj. = 2.90, 2.95 and 3.0, there are three visible 
groups of bins. With increasing, the distribution shifts to 
the right with larger most probable particle number. 



n 


Refs. Uni and [30 More specifically, we consider the path 
integral representation of the thermal density matrix 
(Eq. (6) in Ref. 1331) for the classical XY model. 
However, we use the (ni, 712 , • • • ,nL^) representation and 
the corresponding transfer matrix introduced in Sec. EH 
rather than the original spin variables. In the following, 
the eigenstates of the transfer matrix will be treated as 
quantum states. 

We consider the system, denoted AB, and subdivide 
it into two parts denoted A and B. We first dehne the 
thermal density matrix pab for the whole system 

Pab=T^^/Z. (14) 

We have the usual normalization Tr pab = 1- If the 
largest eigenvalue of the transfer matrix is non degenerate 
with an eigenstate denoted |H), we have the pure state 
limit 

lim PAB = |fl) (HI ■ (15) 

Lt—^oo 

In the following, we will work at finite and will deal 
with the entanglement of thermal state^^. In general, 
the eigenvalue spectrum {pABi} of pab can then be used 
to define the thermal entropy 

St = — ^ PABi (16) 

i 

The subdivision of AB into A and B refers to a 
subdivision of the spatial indices. We define the reduced 
density matrix pA as 

PA = AtbPab- (17) 

We define the entanglement entropy of A with respect 
to B as the von Neumann entropy of this reduced 
density matrix pA- The eigenvalue spectrum {pAi} of 
the reduced density matrix can then be used to calculate 
the entanglement entropy 


FIG. 6. (Color online). Comparison of the particle number 
distribution Pin) from the worm algorithm and TRG with 
different Ds. 


worse by keeping the same truncation dimension Dg. 
When the time dimension Lt becomes large enough, the 
particle density distribution becomes more consistent for 
the cases with different truncation dimension, because it 
is almost centralized in one bin. 


IV. CALCULATION OF THE THERMAL AND 
ENTANGLEMENT ENTROPY 


In this section we explain how to compute the thermal 
entropy and the entanglement entropy in a consistent 
way. We construct the reduced density matrix following 


Se = PAi ln(/9Aj- 


(18) 


The computation of the thermal entropy can be 
performed using the eigenvalues of the transfer matrix 
Xi discussed in the previous section together with 
the normalization pi = Xf*/J2j^f*- Note that 

if we introduce a temperature T and energy levels 
by identifying Af‘ with expi—Ei/T), we recover the 
standard relation St = (E) /T + \nZ. This identification 
makes clear that T cx 1/Lt- Later, in Eq.(19), we will use 
units where T = 1/Lt. 

For the computation of the entanglement entropy , we 
assume that A has 2^^ sites and B has 2^^ sites and 
then perform the Ia and ^B blockings for the subsystems. 
The coarse-graining along the spatial direction ends with 
two sites, one for A and the other for B. We can 
then contract the indices from the two sites in the time 
direction without further truncation. Tracing over the 
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FIG. 7. (Color online). Illustration of the entanglement 
entropy calculation. The horizontal lines represent the traces 
on the space indices. There are Lt of them, the missing ones 
being represented by dots. The two vertical lines represent 
the traces over the blocked time indices in A and B. 

space indices liking A and B, we obtain the transfer 
matrix T^{nA},{n'^},{nB},{n'p}) for the whole AB system. 
Taking the Lf power, tracing over ub, and normalizing, 
we obtain the reduced density matrix { 71 ( 4 }). 

This is illustrated in Fig.jTj 

In the following, we specialize to the case of £ a = £b- 
Assembling subsystems of different sizes is not difficult 
because the space indices are not renormalized in the 
transfer matrix approach, but this will not be done here. 

As a first study, we have considered the cases /3 = 
0.1, Lx = 4 and Lt = 16, 32, 64 and 128 with /r below 
3.2. Both the thermal entropy and entanglement entropy 
develop a peak over the SF phase. We see that for small 
Lt, the thermal entropy is larger than the entanglement 
entropy, but as we increase Lt, the entanglement entropy 
becomes larger than the thermal entropy. The results are 
shown in Fig. For both the thermal entropy and the 
entanglement entropy, a fine structure appears for large 
Lt- This is discussed in the next section. 


V. THE FINE STRUCTURE OF THE SF PHASE 

In this section, we first give an idea of the large 
structure of the phase diagram (where the different SF 
and MI phases are located) and explain in more details 
the fine structure of a single SF phase, more specifically 
along a line of constant (3 and where /r is varied in 
order to interpolate between the the two MI phases with 
successive values of p. 

At small /3, we have mostly MI phases with small SF 
phases in between them. A good qualitative picture can 
be obtained by temporarily considering the anisotropic 
case Px = 0 where the problem is exactly solvable and 
then restoring = p perturbatively. When Px = 0, 

we have only onsite interactions and in the large Lt limit, 
the problem reduces to finding the value n* of n which 
maximizes In{Pi)^^”' for given p and Pp The largest 



FIG. 8. (Golor online). Entanglement entropy (EE, dash 
line) and thermal entropy (TE, solid line) for p = 0.1, Lj, = 4 
and Lt = 16, 32, 64 and 128. 

eigenvalue of T is then {In*{Pi)e^^*)^^. The quantum 
picture is that for large Lt, the relevant state has all Lx 
sites in the n* state and we are in the MI phase with 
p = n*. In summary, the approximate large structure at 
small P is obtained by increasing p from zero and going 
through the MI phases with n*= 0, 1, 2, ... . 

The fine structure of the SF phase between the 
MI phases can be approached by restoring Px = Pt 
perturbatively. The SF phases are approximately located 
near values of p where n* changes. To be specific, 
we will consider the example of P = 0.1, Lx = 4, 
where the transition occurs near p^ = 2.997-•-when 
Px = 0. In this limit, we have 16 degenerate 

states |0,0,0,0), |1,0,0,0), ..., |1,1,1,1) which can 
be organized in “bands” with n = 0 (1 state), n = 1 
(4 states), etc. Below, we call the approximation where 
the indices inside the kets are only 0 or 1 the “two-state 
approximation”. 

The effect of Px is to give these bands a width and lift 
the degeneracy. The energy levels are defined in terms of 
the eigenvalues of the transfer matrix as 

f;, = -ln(A). (19) 

If we plot the energy levels versus p, we see that we have 
successive crossings corresponding to states of increasing 
n. This is illustrated with the two lowest energy levels in 
Fig. [91 Notice the piecewise linear behavior with slopes 
corresponding to the particle number 0, -1, -2, -3 and -4. 
As the levels cross the thermal entropy rises to In 2. 

We observe that near p = 2.90, the lowest energy state 
changes from |0000) to a state with n = 1 

|n,n = 1) = i(|1000)-l-|0100)-f 10010)-f 10001)) . (20) 

It is easy to calculate the reduced density matrix for A 
defined as first two sites and B as the last two sites in 
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FIG. 9. (Color online). The two lowest energy levels (EL) as 
a function of /r for = 4, Lt = 256, /3 = 0.1 and Ds = 101. 
As increases, lines of successive slopes 0, -1, -2, -3, and -4 
are at the lowest level. At each crossing, the thermal entropy 
jumps. The values of thermal entropy are shifted vertically 
by -0.5 to make the figure readable. 


the limit where Lt becomes infinite and for values of /i 
where |r2,n = 1) is the unique ground state 

PA = TrB|0,n=l)(L!,n = l| (21) 

= i(|10)-f |01))((10| + (01|) + i|00)(00| . 

The eigenvalues of pA are 1/2, 1/2, 0, and 0 and the 
entanglement entropy of this reduced density matrix is 
In 2. 

A n = 2 state becomes the ground state near p=2.95. 
It is in good approximation a linear superposition of the 
6 states with two O’s and two I’s. The two states |1010) 
and 10101) have a slightly larger coefficient suggesting 
weak repulsive interactions. In the numerical expression 
of this ground state, we also have contributions from 
states such as |2000) but with a small coefficient. In 
general, for small /3, the two state approximation is good 
(the corrections are small). 

A n = 3 state becomes the ground state near /r=3.03. 
The ground state can in good approximation be described 
as 

\n,n = 3) = i(|0111)-b|1011)-f |1101)-b|1110)) , (22) 


the worm algorithm with Lt=256. Typical world lines 
are displayed in Fig. Given that Lt=256 is relatively 
large, by taking p approximately in the middle of the 
density plateaus (discussed below, see Fig. 111, we obtain 
that most configurations have n corresponding to the 
average value. Figurejl^shows typical results for p =2.93 
(n=l), 3.00 (n=2), 3.07 (n=3) and 3.14 (n=4) using a 
graphical representation similar to Fig. Their spatial 
winding number are discussed in Appendix]^ The most 


/3=0.1;p=2.93 


/3=0.1;p=3.00 




FIG. 10. Worm configurations for p =2.93 (n=l), 3.00 
(n=2), 3.07 (n=3) and 3.14 (n=4) for (3 = 0.1, = 4 and 

Lt = 256. 


which is |n, n = 1) with O’s and I’s interchanged and one 
can interpret the 0 as “holes”. Finally, near p=3.10, 
11111 ) becomes the ground states (with again many small 
corrections). In general, there is approximate mirror 
symmetry about the “half-filling” situation. 

A similar approximate two state description is valid in 
the next SF phase located between the MI phases with p 
1 and 2. One just needs to replace 0 by 1 and 1 by 2. 

We now follow the same path in the phase diagram but 
from the point of view of the word lines generated with 


important feature of these configurations is that almost 
all the |n|’ s are 0 or 1. For a given particle number n, 
between most time slices we have n time links carrying a 
current 1 and — n time links carrying no current. In 
rare occasions, we observe that the line merge or cross. 
Overall, we can think of these configurations as a set of 
weakly interacting loops carrying a current 1. This is 
in line with the dominance of the states like |1010) over 
states like |2000) discussed above. 

We now proceed to calculate the thermal entropy and 
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FIG. 11. (Color online). Fine structure of the entanglement 
entropy (EE, blue), thermal entropy (TE, green) and particle 
density p (red) spanning the MI(p = 0), SF and MI(p = 1) 
phases successively with p increasing at fixed (3 = 0.1. There 
are three different system sizes: i with Lt = 256, = 8 

with Lt = 512, and = 16 with Lt = 1024, respectively. 
The thermal entropy has L^ peaks culminating near ln2 ~ 
0.69; p goes from 0 to 1 in L^ steps and the entanglement 
entropy has an approximate mirror symmetry near half filling 
where it peaks. 


the entanglement entropy for = 4 with Lt = 256 (as 
discussed above), and also for larger lattices = 8 with 
Lt = 512, and L^ = 16 with Lt = 1024. Figure [TT] 
shows the fine structure in the SF phase for increasing 
sizes. In each case, we go through the MI(p = 0), SF, 
and MI(/9 = 1) phases successively as we keep increasing 


p while keeping the other parameters fixed, as already 
illustrated in Fig. We kept the same truncation 

dimension Dg = 101 for the three cases. 

From Fig. we observe an oscillation structure in 
the entanglement entropy and the thermal entropy. The 
number of sites in the spatial direction L^ dictates 
the fine structure. There are L^ transition points in 
the entanglement entropy and Lx peaks in the thermal 
entropy. With Lx increasing, the transition points close 
to the two boundaries become difficult to resolve as 


shown in the case of Lx = 16, L* = 1024 in Fig. 11(c) 


Higher p resolution and larger Ds are needed to obtain a 
clear picture of the oscillations. The approximate mirror 
symmetry of the entanglement entropy with respect to 
the half-filling point persists for larger lattices. As far as 
the thermal entropy is concerned, the height of the peaks 
are all close to In 2, which corresponds to the two-fold 
degeneracy in the ground state. The peaks are located 
at the place where the energy levels from the ground state 
and the first excited state cross as shown in Fig.[^ 

As Lx increases, the fine structure near the two 
boundaries connected to the MI phases is not so 
prominent as the middle regime and a better resolution is 
needed to discern the fine structure. A similar discussion 
applies to the situation when the system experience the 
MI(p = 1), SF, MI(p = j -|- 1) phases successively and 
we have checked that the same type of fine structure 
appears. Note also that for small Lx and Lt large 
enough, p has significant plateaus that could be qualified 
as incompressible regions, however the width of these 
regions shrinks like 1/Lx as Lx increases. 


VI. CONCLUSIONS 

In summary, we have used the TRG method 
combined with the conservation law to calculate the 
particle number density, the thermal entropy and the 
entanglement entropy for the classical XY model with 
a chemical potential. The particle number distributions 
agree well but not perfectly with the results obtained 
with the worm algorithm. The discrepancies seem to 
increase with Lx- This can probably be explained by the 
fact that at small /3, the spectrum has bands with Lj,, 
.Lx(Lx-i) states with increasing n. Some of these 

bands merge when p is increased enough to get in the SF 
phase. As Lx increases, the bands become denser and we 
will typically keep the lowest energy states irrespectively 
of their particle number. Taking into account the particle 
number in the truncation process may help getting more 
accurate distributions. The energy bands are illustrated 
in Fig. 

Besides the large structure associated with the 
alternation between SF and MI phases with integer 
particle number density, we found that in the SF regime, 
there is an interesting fine structure controlled by the 
spatial dimension Lx- The entanglement entropy, the 
thermal entropy and the particle number density vary in 
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FIG. 12. (Color online). Energy levels for /? = 0.1, = 4. 


a way that is consistent with each other. The thermal 
entropy shows peaks located at where the energy 
levels from the ground state and first excited state cross. 
Degenerate perturbation theory explains why the energy 
levels cross while varying the chemical potential. As 
a result, the step-wise structure occurs in the particle 
number density in the SF regime, as already observed 
in Ref. [191 The particle number and the winding 
number of the world lines increase by one unit at a 
time. The entanglement entropy shows steps and an 
approximate mirror symmetry with respect to the half¬ 
filling point. The details of the fine structure depend on 
the ratio L^/Lt and the infinite volume limit of the two- 
dimensional classical model needs to be defined carefully. 

In the approximation where the eigenstates of the 
transfer matrix are made of states with only 0 and 1 
for all sites, or correspondingly if the world lines have 
mostly links carrying currents with |n| equal to 0 or 
1 , this symmetry corresponds to interchanging 0 and 1 
and the particle number n with — n which explains 
the approximate mirror symmetry of the entanglement 
entropy. It is possible that the fermionic picture 
would become more clear if we could establish some 
approximate correspondence with the spin-1/2 quantum 
XY model. 

We also would like to mention some analogies. As the 
chemical potential can be interpreted as an imaginary 
gauge field in the time dir ection , it is not surprising that 
studies of Polyakov’s loopP^^^ (Wilson loops closing in 
the periodic time direction) show similar fine structure. 
Note also that the crossing pattern of the energy levels as 
a function of ^ found here resembles the energy crossing 
found in a study of the spectrum of rotating tubes as a 
function of the rotation rate for a proposed cold atom 
experiments^. 


Appendix A: Remarks about Fig. 

In this appendix, we discuss the sign convention and 
the spatial winding number of Fig. We also explain 
why this configuration is typical. 

The sign of the n’s associated with links with |n|=l can 
be figured out from the following information. All the 
(vertical) temporal link indices are positive. The sign 
of the (horizontal) spatial links can be obtained using 
the conservation law. On time slices 5, 7 and 17, the 
current moves to the right and the sign is (by convention) 
positive. On time slices 21 and 22, the current moves to 
the left and the sign is negative. In this configuration, 
the winding number in the temporal direction is 2 and 
the winding number in the spatial direction is 0 (there 
are as many positive as negative spatial n’s). 

The fact that the configuration of Fig. [^is typical for 
/3=0.1 and /r=3 can be understood from the numerical 
values of the weights. Because of the large value for fi, the 
weight for the temporal links with n=0 (/o(O.l) ~ 1.0025) 
and n=l (/i(O.l) exp(3) ~ 1.00553) are almost the 
same, while the weight for n=2 (/2(0.1) exp(6) ~ 0.5047) 
is smaller. The weight for n=-l (/i(O.l) exp(—3) ~ 
0.002492) is very small and there are no temporal links 
with negative values of n in the configuration. For 
the spatial links, the relative cost of a lateral move is 
/i(0.1)//o(0.1) ~ 0.05 and there are only 6 lateral moves. 
The fact that there are only 2 temporal links with n=2 
can be understood from the fact that the merging of two 
n=l lines into one n=2 line requires one lateral move in 
addition of a weight about twice smaller. 


Appendix B: Spatial winding numbers in Fig. |10| 

In this appendix, we discuss the spatial winding 
numbers of the configurations of Fig. For /r=2.93, 
all the time (vertical) links have n=l. For the spatial 
links there are 20 right movers and 16 left movers (so the 
spatial winding number is I). For /r=3.00 (n=2), between 
most time slices, there are two vertical links carrying a 
n=l current, and the two lines only merge 4 times into a 
single n=2 line for a small number of time steps (the total 
is 5 vertical links with n=2). There are 26 right movers 
and 34 left movers (so the spatial winding number is -2). 
For II—3.07 (n=3), between most time slices, there are 
three vertical links carrying a n=l current. There are 5 
occurrences where two n=l merge into a single n=2 line 
for a small number of time steps (the total is 10 vertical 
links with n=2). There are also 3 crossing (points where 
the four lines attached all have |n|=l). There are 19 
right movers and 15 left movers (so the spatial winding 
number is -1-1). For ^=3.14 (n=4), 13 = 0.1, = 4 

and Lt = 256, there are essentially four parallel vertical 
lines each carrying a n=l current except for 4 occurrences 
where they briefly merge (11 n=2 vertical lines, 4 right 
movers, 4 left movers, no spatial winding number). 
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